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APPLICATION OF INTEGRATION ALGORITHMS IN A PARALLEL PROCESSING 
ENVIRONMENT FOR THE SIMULATION OF JET ENGINES 


Susan M. Krosel and Edward J. Milner 
National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 


Abstract . The development of digital dynamic simulations requires careful selec- 
tion or an appropriate integration algorithm. This paper illustrates the appli- 
cation of predictor-corrector integration algorithms developed for the digital 
parallel processing environment. The algorithms are implementca and evaluated 
through the use of a software simulator which provides an approximate representa- 
tion of the parallel processing hardware. Test cases which focus on the use of 
the algorithms are presented and a specific application using a linear model of a 
turhofan engine is considered. Results arc presented showing the effects of 
integration step site and the number ol processors on simulation accuracy. Real- 
time performance, inter-processor communication, and algorithm startup are also 
d i scussed . 


INTRODUC TION 

In recent years, there ha? been a growing need to obtain simulation models 
which represent physical systems over tb*;r entire operating range. An example 
of this is the jet engine. Th" > demana for higher performance in these system.", 
has resulted in increased system complexity and a need for more in-depth analysis 
of their dynamic behavior. There is an additional need for detailed system mod- 
els to support the development of digital controls for these systems. In both 
the design and evaluation of these controls, simulations are frequently used. 

Digital computers are used extensively for simulation because of their ease 
of programming, repeat ibi 1 i ty of results, and to a large degree, the portability 
of the simulations. Digital s imul at *. ons , such as GENENG, DYNGEN, and NNEP (Ref- 
erences 1,2,3) provide the capability of predicting the steady-state and dynamic 
performance of a wide variety of gas turbine engine configurations. Digital com 
puters, however, are limited in their usefulness for time-critical simulation 
applications by their inherent sequential execution of program instructions and 
serial computation within these instructions. In applications such as the vali- 
dation of digital control hardware and software, the requirement for real-time 
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response of the simulation has necessitated either the uae of large, dedicated 
computer systems with instruction cycle times in nanoseconds or the simplifies* 
t ion of the model . 

With the advent of and current advances in digital micro-computer technology, 
it is now possible to develop small, compact, computer systems for simulation. 

More importantly, it may r..iw be possible to implement a detailed simulation model 
and still achieve real-time operation. This will permit the aimulatlon to be 
used in a wide variety of applications including digital control system develop- 
ment, checkout, and troubleshooting as well as performance studies. One approach 
that has been proposed is to connect several microprocessors in a parallel ar- 
rangement and to provide a means of communieat ion between the processors. The 
simulation is then partitioned over the several processors by dividing the syatem 
equations among the N p ocessors forming the parallel digital system. However, 
partitioning necessitates a careful and thorough consideration of the dynamic 
coup 1 within the model to determine the optimal breakdown of the system func- 

tions. In some cases, inherent parallelism in the system may simplify the parti- 
tioning. The issue of how many processors to use then can be addressed. For 
efficient operation, the portions of the simulation that are allocated to the 
individual processors should use approximately the same amount of compute time 
per processor. This will insure crrrect updating of system variables and avoid 
wasted time in the calrulet ion cycle. The updating of variables within the par- 
titioned simulation will require not only careful timing considerations but also 
efficient data transfer between processors to avoid inadvertent phase shift. 

The development of digital simulations, in general, depends on the selection 
and implementation of suitable numerical integration algorithms. These algor- 
ithms should provide for accurate and efficient solution of the differential 
equations that describe the svstem being simulated. For the gas turbine engine, 
these equations are typically nonlinear and involve multivariable function* that 
describe the performance of the engine's rotating components (fan, compressor, 
and turbines). In general, most of the computing time is used in the calculation 
of the system derivative function. Therefore, in the selection of an integration 
algorithm in a time-critical application, one must consider the number of deriv- 
stive calculations associated with the algorithm. Much work has been done in the 
study of integration methods for e single processor system (Reference* 4,5). The 
design and application of integration algorithms for a parallel-processing system 
depends on additional factors, such as: the number of processors, the method of 

partitioning the simulation, the inter-computer data transfer mechanism, the 
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c 'inputs t ional speed of the processors, end the need to input end output simulation 
t ta (References 6,7). For example, it may be possible to partition a problem 
into linear and nonlinear parts and to use different Integration algorithms on 
each part (Reference 8). 

This paper discusses the application of parallel predictor-corrector 
algorithms (Reference 7) to the simulation of a typical turbofan engine. The 
simulation is Intended to run on an MIMD (multiple instruction > multiple data) 
parallel processing system (Reference 9). A software simulator was used to pro- 
vide an approximate representation of a parallel processing system and to evalu- 
ate the performance of the algorithms. 

A first-order system and a second-order system were used to evaluate the 
algorithms. Each of these systems was excited by a unit step. The effects of 
the number of processors used and the integration stepsize on simulation accur- 
acy, resolution, and stability were determined. Results are presented and dis- 
cussed in the following sections. 

AN ALTERNATIVE TO PARTITIONING 

Miranker and Liniger (Reference 7) suggest a parallel predictor-corrector 
integration technique which merits consideration as an alternative to partition- 
ing. No partitioning i6 required because the algorithms they present require 
that the entire simulation reside on each processor. Normal!;* when using 
predictor-corrector integration, the current corrected value of a parameter is 
based on its current predicted value. The operations are sequentially regimented 
with the requirement that calculation of the current predicted value be completed 
before calculation of the current corrected value may begin. However the 
Miranker and Liniger algorithms predict and correct current values based on val- 
ues already evaluated in a previous calculation cycle. Hence, prediction of some 
values and correction of others can take place simultaneously. 

This charcter 1st ic of the algorithm allows concurrent calculation to take 
place on parallel processors. Taking advantage of this assumed calculation 
power, the algorithms are able to eperate on more than one integration time step 
during a single computer simulation update cycle. The update cycle time remains 
fixed; but since more than one Integration time step, h, is calculated during 
this period, the simulation is effectively speeded up. Ideally, unless stability 
or accuracy problems arise, using a sufficiently large number of parallel proces- 
sors should allow real-time operation. The relationships among the number of 
processors, system stability, and system accuracy />rc- examined in detail in the 
RESULTS AND DISCUSSION section. 



Consider a dynamic system described by equations of the form 


y* • ffx.y), k > 0, yfO) - y 0 <1> 

Suppose 'bat numerical solution is required at the mesh points * n ■ in - l)h, 

n - 1,2 where h is the step size. Hiranker and 

this system by using predictor-corrector algorithms such 
order, four-processor algorithm 

y 2n+2 " y 2n-2 4 * hf 2n 
y 2n+l " y 2n-2 4 ? ( f 2n 4 f 2n-l) 
y 2n ' y L * 2 ( 3f 2n ' 9f 2n-l) 
y 2n-l ” y 2n-3 4 2hf 2n-2 

where v is an approximation to y at x , y^ is the 

■ n r n n 

x^; Yp is the corrected value at x n , » f^x n ,y^ and 

Pred ictor-correetor algorithms for one, two, and four processors were included in 
Reference 7. An eight-processor algorithm was derived by the authors based on 
the Miranker and Liniger methods. Second-order algorithms used for this study 
ere listed in lable 1. 

Figure 1 shows the allocation of the parallel predictor-corrector algorithm 
equations for the four-processor case. As shown in Figure 1, the algorithm for 
four processors calculates predicted values on two of the processors (A and b) 
and corrected values on the other two (C and D) . Thi6 results in two outputB 
(i.e., corrected values) per calculation cycle. The outputs for the four pro- 
cessor algorithm are, in order, the corrected values from processor D and then 
from processor C. During this same calculation cycle, the predicted values are 
being calculated on processors B and A for one and two output times in the 
future respectively. These values then feed back into the algorithm for the 
updating of the corrected values in the next calculation cycle. In each new cal- 
culation cycle, the current values of the inputs are brought into the calcula- 
tions on each processor. For each state of the system being simulated, eleven 
transfers of data mu6t be accomplished for the algorithm for each step In the 
integration cycle. It should be noted that only N- 1 of the derivative function 
calculations are actually used in the N-processot algorithm. 


Liniger propose solving 
sa the following second - 
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(3) 

(4) 

(5) 
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SIMULATION APPROACH 

To test the predictor-corrector algorithms before Implementing them on 
actuel parallel processor hardware, a software simulator was developed. The 
software simulator ; s written in Fortran to run on an IBM TSS 370. The program 
is structured in a modular fashion through the use of subroutines (t igur-* 2). 
Subroutines are used to represent the N processing elements. Psuedo- 
parallelism is achieved through the use of argument lists for variable transfer 
with distinct variable names during a calculation cycle. The code for the cal- 
culation of the system derivatives is contained in one subroutine. This permits 
enRy modification of the software simulator for different systems being simu- 
lated. The software simulator accurately represents the problem-solving phase oi 
the parallel processing predictor-corrector integration algorithms. Careful 
attention has been given to the outputting ol results with respect to time. 
Representation of actual simulator control was not incorporated into the software 
simulator. Coding to represent a simulator controller and an input-output pro- 
cessor has been included at a very simplistic level. 

RESULTS A M 1 DISCUSSION 

The Miranker and Liniger integration technique was examined from the real- 
time simulation point oi view by applying it to a linear first-order system, twe 
diiferent linear second-order systems, and a linear fourth-oruer engine model. 
Following is a detailed examination of each. 


M RSI -ORDER SYSTEM 

The response of a first-order Jag to a unit step input obeys the relationship 

iy' + y - 1 (6) 

The closed-form solution to Equation (6) with the initial condition y ( 0 ) * 0 is 

>’ ■ 1 - exp(-x/T> (7) 

This response which is dependent on the value of the system time constant, t, 
can be displayed as a single curve provided that the parameter x/t is con- 
sidered the independent variable. Rearranging Equation (6) yields 


Thus , 


y’ - { (i - y) 


• f ( t.y) - { g(y) 


<B) 


y* 


(*> 



Hcncr, the second-order four processor predictor-corrector algorithm can b« 

written of 


■ v 2n4 2 * y 2n-2 4 4 7 *2n 

(10) 

y 2n+l “ y 2n-2 4 7 t (*2n 4 *2n-l) 

(11) 

c „ V C 1 h L P Q P \ 

y 2n y 2n-3 ' 7 7 \ 3 *2n ‘ 9 2n-lJ 

(12) 

v C - v C + 2 h e C 

y 2n- 1 y 2n-3 4 1 t *2n-2 

(13) 


It is clear that the predictor-corrector solution (sequence of points) is depend- 
ent nn the number of processors, the system time constant, t, and the integra- 
tion step size, h. For a fixed number of processors, the result is a family of 
solutions each corresponding to a different value of the parameter, h/t. 

To determine how well the predictor-corrector solutions jeompare with the 
closer-form solution (Equati on (7)), a figure of merit was established in the 
following manner. The absolute value of the difference between the predictor- 
corrector solution ana the closed-form solution was integrated over the time 
interval 0 to iti. This cumulative error was then expressed as a percentage of 
the total area bounded by the closed-form solution over that same time interval. 
This can bi rypresseu by 



(closed form - algorithm|dx 

r 


r ' 

J Jcloseo forn((!x 


(14) 


Thus, each of the predictor corrector solutions has associated with it a relative 
error . 

In studying the effects of varying the number of processors and/or the inte- 
gration stepsize, it is helpful to define the "effective step advance per calcu- 
lation cycle," H, as follows 

H - h*(N/2) (15) 

where N is the number of parallel processors used. If only one processor is 
used, two derivative calculations are required per integration time step. One 
derivative calculation is required for the predictor and another one is required 
for the corrector. Consequently, the effective computation time for this case is 
twice the calculation time and the effective step advance per calculation cycle 
is b/2 . in the four-processor case, only one derivative calculation per pro- 
cessor is required and each calculation cycle will result in two time steps being 
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t- ;• 1 1 i'l,n t>(> . Thus the effective step advance per calculation cycle 1# 2h. In 

Muuvinr tin’ effect of changing the number of processors or using different inte- 
gration step fi7>s, comparisons should be made on the basis of the same effective 
f t cp advance per calculation cycle (that is, on the basis of equal values of H). 

It was previously noted that the parameter b/t could be used to elimi- 
nate the effect of the system time constant, t, on the predictor-corrector 
solutions for a fixed number of processors. Similarly, the use of H/t allows 
results from varying numbers of processors to be compared. The percent error 
resulting from using the second-order predictor-corrector algorithms as a func- 


tion of the paremeter H/t is shown in Figure 3. Data corresponding to one, 
two, four, and eight processors are presented. We see ‘.hat, for small values of 
H/t, there is little, if any, advantage to using multiprocessors. A Single 
processor gives good accuracy and does not require the transfer of any data. 
Hcncr , the mechanics of the simulation is kept simple. As the value of H/t 


increases, v.e sec that using a larger number of processors gives better accu- 
racy. Fot a value of H/t • 0.121». accuracy is improved by almost a factor of 


four by using eight processors. This accuracy improvement is due in part to the 
larger number of points obtained from using more processors (a finer grid). When 
using eight processors, four points are calculated per calculation cycle. The 
effective step advance using eight processors consists of three intermediate 
points plus an end point. When using only two processors, one point is calcu- 
lated per calculation cycle. In essence, then, when using two processor , all wc 


arc calculating in our solution are the end points of the eight processor ca'e; 
and, hcncc, we are suffering a loss in accuracy by using a smaller number of 


processors. 

As H/t approaches a value of 0.25, we see that slightly improved a r cu- 
raev is obtained using four processors rather than eight. This may be due to the 
step s i 7 e becoming too large for the more complicated eight processor algorithm. 
For H/t equal to 0.40, the eight processor algorithm is unstable. The inte- 
gration step size has become too large for this complicated algorithm to hold 
together. However, the four processor algorithm is still stable at this point. 
This brings out an important point; namely, a simulation is not necessarily 
improved by using more processors. 

For the purpose of illustration, suppose that a real-time simulation of a 
first-order system is desired. Suppose further that the corner frequency of the 
system is 50 radians/sec, that four processors are available for the simulation. 
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and that each processor takes 2 milliseconds to compute the required derivative 
function. 

Hence, calculation time, t £ , is 0.002 sec and the time constant, 
t * 'corner frequency)* 1 , is 0.02 sec. Since t>ol-tim« simulation is desired, 
the effective step advance must be equal to the calculation time, t c . Thus 
H » t c or H/t ■ t c /t ■ 0.002/0.02 - (.10. From Figure 3 we see that the corre- 
sponding simulation erroT will be approximately 1,3 percent if all four proces- 
sors are used. Let resolution be defined as the calculated number of points per 
cycle at the highest frequency of interest. Then, resolution for this case is 
approximately 125 points/cycle at 50 radians/sec. Also from Figure 3 we see that 
the simulation error will be approximately 3 percent if one processor is used. 
Resolution for this esse is still sufficient wi.h approximately 30 points/cycle 
at 50 radians/sec. Now suppose that calculation time t c increases from 2 milli- 
seconds to 3.6 milliseconds. If 2 percent is the maximum allowable error and if 
four processors are available, we see from Figure 3 that a maximum value of 
H/t - 0.13 is permitted. However, for t c - 0.00.6 sec and t • 0.02 sec, 

H/t - 0.16 for real-time operation. Hence, the simulation cannot run real-time 
under these conditions. 

Whac can be done to correct the situation? There are several alternatives. 

If eight processors were available, the simulation would run in real-time with 
2.1 percent error with H/t - 0.18. If more than four processors were ?ot 
available, then faster processors would be needed. If the hardware was fixed, 
some of the higher frequency dynamics of the simulation would have to be sacri- 
ficed for the aake of real-time operation. A time constant of t ■ 0.028 sec, 
which corresponds to a corner frequency of 36 radians/sec, would be required for 
real-time simulation with 2 percent error using four processors (for then 
H/t • 0.0036/0.28 * 0.13). In determining the numbe- of processors to use to 
simulate a system with first-order characteristics, then, the choice depends on 
several factors. As we have seen, going to more processors gives better resolu- 
tion end can give better accuracy. However, going to more processors may de- 
stabilize an otherwise stable simulation. 

SECOND- ORDER SYSTEM 

As a further teat of the predictor-corrector s'gorithms, a tecond-oroer 
system consisting of a first-order lag feeding another first-order lag was simu- 
lated. The transfer function for such a system is 1 /((tjS + 1)(tjS + 1)) and 
its associated differential equation has the form 


V 


v ? y " 4 (t i T T 2 )y ' * y m 1 

Fit the initial conditions v(0) * 0 ano y'(0) • 0. Equation (lb) hat tha 
c 1 o*eri- f orm solution 


y * 


V "* /T i - 


’i 


t 2 e* x/, 2 

~2 


U6) 
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For the study, two different cases were considered CASE I has the values 
» 0.01 sec and t 2 * 0.02 sec; and CASE 11 ha6 the values ij ■ 0.002 sec 
and ij, » 0.02 sec; the procedure followed was similar to that used for the first- 
order sv6tetn to obtain date. The second-order system responses can be reduced to 
a family of curves by using x/tj as the indepenoent variable. In this case, 
however, each curve represents a particular value of the ratio T 2^ T 1‘ ^ or 
two costs, <rp/tj - 2 and tp/’j “ 10) solutions for varying integration step 
size, h, were obtained using the predictor-corrector algorithms for one, two, 
four, ami eight processors. As was done. for the first-order system, an error 
measure was calculated for each solution. The time interval consioered for this 
computation was 0 to four times the smaller time constant, tj. The error data 
for these cases were presented ns a function of H/rj to facilitate comparisons 
of results obtained fmm different numbers of processors. 

Figure U presents the percent error versus H/ij for CASE I. This system 
h..s the slower dynamics cf the two cases considered * 100 racians/sec ano 

»»o • 50 rndisns/sec) . Notice that i he er r or obtained using either one processor 
or two is virtually identical. However, using two processors gives better resolu- 
tion (at least 15 points per cycle at the highest frequency compared with 8 points 
per cycle using one processor). Notice also that, for a value of H/ r j • 0.5, 
one processor is stable and two or more are not. 

If four processors can be used, improved accuracy and improveo resolution 
are realized. At least 30 points per cycle arc obtained at the higheat fre- 
quency. Notice also that, at H/tj * 0.t(, accuracy is improved 100 percent. 

Figure U shows that there is no reason to use eight processors. There is only a 
narrow range of stability ior eight processors anc both accuracy arc, resolution 
over this range are excellent with four processors. 

Percent error versus H/t] for CASE 11 is presented in Figure 5. This 
system has faster dynamics ■ 500 radians/sec and * 2 * 50 radians/sec) than 
CASE 1. CASE 11 alio has a much smaller range of H over which the simulation 
is stable due to the higher frequencies involved. (The time constant ij is 
five times, smaller than in CASE 1.) However, accuracy is excellent over the 
stable range, no matter how many processors are used. As in CASE 1. the error 
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« M aincd using either one processor or two is virtually identical? howavar battar 
msolut ion obtained with two proeeRSots due to the finer spacing of points in 
the solution sequence. Resolution at 500 radians/aec, using one processor, is 
approximately b points/cyclr • using two processors, it is approximately 12 
points/cycle . Because of the excellent accuracy and resolution over the stable 
range for this caae, there should be no need to use more than two processors. 

Attaining real-time operation for a stiff system such as this would b« 
difficult. The »j high frequency is a dominant term requiring a very small 
integration time step for stability of the simulation. Since for real-time 
simulation t c * H, we see that this requires the calculation cycle time to be 
;>o more than a millisecond. This calculation time may be met for elementary 
systems with available computers. However, »>r simulations of complex systems 
such as turbofan engines, to achieve real-time operation probably will require 
sacrificing some of the higher frequency dynamics. The most favorable condition 
for real-time simulation in the case is with one or two processors. 

ENGINE MODE). 

As a final test to determine whether the parallel predictor-corrector inte- 
gration a\f r i thms are applicable to the turboian engine simulation problem, a 
state-space model of a representative ei.gine was used as the system in the soft- 
ware simulator. 

The model selected was a reduced-order (fourth-order) linear model at the 
sea-levr], static, intermediate power operating condition. It was obtained from, 
a full state (Ibth-oraer) linear model by normal reduction techniques (Reference 
10). This linear model was validated along with other linear mocels of the 
engine in Reference 10. The reduced-order model still retains important dynamic 
character i st ica of the system but is easier to handle mathematically. The 
'mat hemat i cal representation o ( the system is given by the following equation* 


A x 

♦ B u 

(18) 

C 5? 

4 15 U 

(IV) 


where Equation (18) is a linear, constant-coei f icient-metr ix differential equa- 
tion (valid only at a given operating condition) that represents the computation 
of state variable derivatives. Matrices A, the system matrix, and B, the con- 
trol matrix, show the sensitivity of the time derivatives oi the state variables 
T to variations in the state variables T and control inputs IT. Equation (IS 1 ) 
is a linear, conatant-coeff icient matrix algebraic equation that represent* the 
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eomputat ion of observed engine parameters. Matrices C, the output matrix, and 
D, the d irect- couple matrix, relate the changes in the observed perimeters ~ to 
the variations in th? state variables x and the control inputs u. For the 
selectee reduced-order mooel, the states, represented by vector T, are fat' 
speed, compressor speed, compressor discharge pressure, and interturbine pres- 
sure. The observed parameters, represented by vector J, are engine net thrust, 
tr’al engine airflow, hutner-exit temperature, fan stall margin, compressor stall 
margin, measured fan-exit 4p/p parameter, and calculated fan-exit sp/p 
parameter. The control inputs, represented by vector IT, are main burner fuel 
flow, exhaust nozzle Jet ares, fan inlet guide vane position, compressor variable 
vane position, and compressor bleed flow fraction. 

The mortices A. B, C, ano D are giver, in Table ,1. ho single linear 
nuclei can accurately represent an engine over its entire operating range. There- 
fore, many linear tnooels art’ tvpically derived at various flight conditions and 
power settings throughout tiic engine operatirg envelope. These models would be 
connected together in some manner tc form a more accurate and representative sim- 
ulation oi the engine process. The selected engine model was derived at the 
sea-level, static, intermed iate power operating condition enu hence is valio only 
in that region. 

To determine if the parallel predictor-corrector alogrithms are suitable and 
appropriate integration methods for the engine simulation problem, a 3-percent 
step in fuel flow was input to the reduced-order model for one, two, four, ano 
eight processor pred i ct ir-correc f or algorithms. The resulting transients were 
compared with an exact solution obtained through evaluation of the state 
t rans it ion matr ix . 

In running the transient on the simulated multiprocessor systems, it was 
found that, as the number oi processors increased, the integration step size for 
the algorithm had to be decreased. This decrease reflects the loss in stability 
due t-> these algorithms as more proceaaors are usee. Table 111 gives the tabula- 
tion of the number of proceaaors and the maximum allowable timestep for etebil- 
ity. These results are summarized in Figure 6 which showa the eilect of the num- 
ber of proceaaora on the required step aize for atable operation. Data are pre- 
sented for the firat-order system, the second-order systems, and the engine 
model. In the engine model, Table II aiowi that the system possesses widely 
spaced eigenvalues, requiring the use cl small tjmtsteps for stability at the 
high freouency dynamica. However, Table 111 ahowa that the algorithms themselves 
require the use of smaller tjmesteps to asrure the stability of the algorithm 
If real-time operation ia desired, the decrease in the integration timestep means 
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t V. ,• 1 t h» Mp,i a 1 low.iblr to compute that system derivative. also decreases. Table 
111 f]*r gives the time allowed for the derivative calculation for real-tie* 

operat i( - . 

The transient response of two atatea •• fan apeed and compressor discharge 
trrssure - and two observed parameter* * burner exit temperature and compressor 
stall marg1r - are shown In Figure 7 for the one, two, four, and eight processor 
ripe? err the erect polutior. The transient 1* shown only for the firat aecond 
apc. it ip aeon that agreement it quite gooo . After the firat half aecond, the 
difference between the exact and any of tne Crtea become* In* igni f leant . The 
on 1 v noticeable error ia aeen in the firat tenth aecond of ';h* transient. Thia 
is ..hewn mr re clearly in Figure 8. Thia error at the beginning of the transient 
i* caused by the startup delay inherent in the predictor-corrector algorithm. In 
»rnrr,-] pr ; <' i r * nr -correct or integration algorithms are quite accurate provided " 

M wbl e timestep it chosen. This ha* been shown in Figure 7. However, the cis- 

i 

linct r isadvartagr that these algorithms possess is in their non-self -start ing 
feoturr. This was shown in Figure 8 by emphasizing the firat 0.1 second of the 
transient rrsponee. This is eounterbal anceo by the need for fewer evaluation* of 
t be svstem derivative, thereby requiring less computation time for the algor- 
ithm overall. If non- real - 1 ime operrt'en is acceptable, then these parallel pre- 
i: ic t or- ror n c t or algo-ithms wav be useful. 

Figure S gives a series ol eight processor transientr tor varying timesteps 
f( 0.00b'. (J . 0002 !> , 0.0007, and 0.0001 aecond over a 0.5 s*-ond range. It can be 
seen trat ar. integration step sire (i 0.0005 second produces oscillation* in the 
(irst O.i second of the transient which subsequently *ii» out. This it due to the 
instability ol the alogrithm at that step. The eight processor algorithm ia 
stable when h is reduced to 0.00025 second. This implies that for real-time 
operat ior with this model, the evaluation of t Vie derivative Junction muat be com- 
pleter, in 1 m.llisecond (see Table 111). Use of a smaller timeatep doe* reduce 
these oscillations but results in tJmesteps much too small to have practical ap- 
p 1 i ca t ion . 

Figure 10 ahows a comparison of the predictor-corrector and actual response- 
over the first tenth of a second. 


'•phOLUSlO t 5 

The Wiranker and Liniget predictor-corrector alogrithma can eliminate many 
problems associated with digital multiprocessing. The meat obvious eevantage i* 
that the predictor-corrector alogirhtms do not require partitioning of the simu- 
lation mocel. Jn applying the algorithms to a linear first-order system and t* 
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two different linear record-order systems, the errors tbtained using one proces- 
sor and two processors were virtually identical. In addition, using four proces- 
sors generally cut this error in half. The range of stable timesteps was de- 
creased by using more processors. In fact, in some cases, increasing, the number 
of processors had a destabilizing effect. Using more processors did give better 
resolution but, in many cases, at the expense of decreasing the stability of the 
simulat ion. 

Use of parallel predictor-corrector integration algorithms poses some dis- 
advantagp when real-time operation is desired. It has been shown that the u6e of 
these algorithms with increasing numbers of processors does result in a reduction 
in the stability of the algorithm and requires the use of smaller and smaller 
integration timesteps. This requirement of very small timesteps means that the 
computing device may have to perform more calculations over a specified time in- 
terval and be able to calculate the system derivative function in a shorter 
amount of time. 

A1 so it has been seen that che predictor-corrector algorithms in themselves 
are not sel f - star t ing , therefore, either a deedt ime occurs at the beginning of 
the transient (no response from the system) or some means of starting the algo- 
rithm must be devisee. Since the calculation of the derivative is only required 
cn N-l processors (ior the N-processor algorithm), this ‘free calculation time’ 
is available for implementation of a starting method. Proper switching logic 
must be incorporated to insure that the predictor-corrector algorithm tales over 
at the end of the startup period. 

The parallel predictor-corrector integration algorithms may be a possible 
means to attain real-time operation for dynamic system simulation if: 

(1) Microprocessor internal cycle times continue to decrease which will 
allow for faster calculation of the derivative functions; 

(2) The cost of microprocessor memory continues to decrease allowing for 
the acauisition of enough memory for algorithm implementation; 

(3) The system to be simulated does not have widely spaced time constants 
or very high dynamic frequencies such that the stable range of time- 
steps is not critically limited; 

(d) Simplification of the simulation model, when necessary, can be made 
without excessive less of accuracy. 

In situations where parallel predictor-corrector integration algorithms can 
not be used, other approaches such as problem partitioning and implicit integra- 
tion shouia be considered. 
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TABLE I . - PARALLEL PREDICTOR-CORRECTOR INTEGRATION 
ALGORITWS (ORDER 1VJD) 


One processor: 


Vl 


* 7 ( 3, n • 4l) 




TVo processor : 


4) - 4l * 2 < 




Four processor: 

y 2n+2 “ y 2r 2 4 


+ 4hfo_ 


P C . 3h / f P , f P \ 

v 2n4] ' y 2n -2 4 T \ 2n 4 f 2n-l/ 


4 ■ v 2n-3 - 7 ( 3 4 ' ,f 2n-l) 


4-1 ‘ 4-3 * 2b 4-2 


Eight processor: 
P 


v, , c # + 8hff 

- 4r>+4 4n-4 4n 


P C 7h / r P , f P \ 

y 4m3 " y 4n-4 + T \ f 4n 4 f 4n-l / 


v 4r>+2 " y 4n-4 4 fthf 4n-) 


P C 5h / f P , f P \ 

y 4rr+] y 4n-4 Y~ \ 4n 4n-l/ 


C. C 5h P c f P 
y 4n y 4n-5 \ 4n 4n-l ( 


( 3 4 - 5t Ll) 
( 3 Ci - st L 2 ) 

r ( 3f <m -2 • 5f 4i>3) 

in- 3 " 4-7 * ‘"(’t -3 ' 2 4-I.) 


C c 5h 
y 4n-l " y 4n-6 ‘ T 


C C 

y 4n-2 * y 4r»-7 
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TABLE 11. - 

ENGINE MODEL (4 f h ORDER, SEA- LEVEL- STATIC, 1NTE»€D1ATE) 

Model equations: 

X 

-Ix+lZ 





y 

* C * ♦ B u 



Mtxiel matrices: 







System matrix, K 




1 

2 

3 

4 


1 

-3.684 

-0.6320 

103.0 

-789.0 


2 

1.200 

-6.137 

79.95 

-522.8 


3 

-1.659 

5.819 

-154.8 

729.0 


A 

0.2584 

-0.2028 

5.581 

-103.7 




Control matrix, t 




1 

2 

3 

4 

5 

1 

0.7260 

-257.7 

-107.9 

-1.972 

-0.1500D 05 

2 

0.6320 

-458.4 

13.40 

-73.93 

-0.13270 05 

3 

0.8555 

647.3 

-10.04 

59.29 

0.7941D 05 

4 

0. 1891D-01 

-1229 

10.16 

-1.938 

-399.6 



Output matrix, C 




1 

2 

3 

4 


1 

1.109 

-0.8876 

24.65 

70.15 


? 

0. 1382D-01 

0.3) 52D-05 

-0.7751D-07 

-0.9650D-02 


3 

0.1093 

-0.1181 

3.136 

-24.47 


A 

0.7315D-0A 

0.5169D-05 

0.5422D-06 

-0. 1508D-01 


5 

-0 . 3066D-04 

0.1294D-03 

-0.2594P-02 

0.1449D-01 


6 

0.4003D-04 

-0.3016D-04 

0.1166D-03 

-0. 1430D-01 


7 

0.1332D-04 

0.4075D-05 

-0.9470P-05 

-0.3745D-02 




Pircct- 

■couple matrix, 

D 



1 


3 

4 

5 

1 

0.1647 

-232.1 

41.22 

-8.291 

-3459. 

2 

0.6344 P-04 

0.1782 

0.7403 

0.2073D-03 

0 . 7989D-02 

3 

0.1329 

-23.97 

1.113 

-1.199 

-2464. 

A 

0.3239D-05 

-0.3606D-0? 

-0.3220D-02 

0. 67 36 D- 04 

0.1239D-01 

5 

-0.8887D-08 

0.1247D-01 

-0.29150-03 

0. 2437D-02 

-0. 1710D-01 

6 

0.331 ID-07 

-0.1 163D-01 

0. 17370-02 

-0.3235D-03 

0.5144D-02 

7 

0.4271D-07 

-0.2997D-0? 

0.7025D-03 

0.5175D-04 

0.2557D-02 


where: 

The states represented hv vector T, arc 
x] Fan speed 

x2 Compressor speed 

x3 Carpresscr discharge pressure 

x4 Tnterturbine pressure 

The output variables, represented by vector y, are 
yl Bigine net thrust 

y2 Total engine airflow 

y3 Burner-exit twperature 

y4 Fan stall margin 

y5 Compressor stall margin 

y6 Qnpirical fan-exit Ap/p parameter 

y7 Theoretical fan-exit Ap/p parameter 
The control inputs, represented oy vector u, are 
ul Main-burner fuel flow 

u2 Exhaust-nozzle- jet area 

u3 Fan- inlet-guide- vane position 

u4 Compressor variable-vane position 

u5 Compressor bleed flow fraction 

Model eigenvalues (sec'^) Corresponding time constants (sec) 


Real 

Imaginary 


-3.1516863 

0.000000 

0.3173 

-5.5467604 

0.000000 

0.1803 

-60.662508 

0.00000 

0.0165 

-198.91636 

0.00000 

0.0050 



1 ABU 111. - CORRELATION AMONG NUMB* OF PROCESSORS, MAXIMUM 1N1EGRATI0N 
STEPS1ZE FOR STABILITY. SYSTEM DERIVATIVE CALCULATIONS. OUIW1S PER 
CYCLE , AND ALLOWABLE TIME FOR DERIVATIVE CALCULATION FOR REAL-TIfC 


— 

Ntmbcr ot 

processor* 

N 

Maxiimrr 
Integration 
stepsize for 
stability 

Nnex 

aec 

Nunber of 
ays ton 
derivative 
calculations 
per processor 
D 

Number of 
outputs 
(corrector 
value*) 
per cycle 
C 

Allowable tine for 
ayaten derivative 
calculation for 
real-tine operation 
Riven by N/2 + ^ 

T . 

aec 

1 

0.005 

2 

I 

0.0025 

2 

.0025 

1 

1 

.0025 

u 

.001 

1 

2 

.002 

8 

.00025 

1 

4 

.001 



PREDICTOR 

EQUATIONS 


CORRECTOR 

EQUATIONS 


Figure L - mvtonantatlan of four praceisor parNM predictor -corrector (ntogrtokm (tgorlthm. 
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Figure 4 - Relative error lor second order system 
1/lTjS ♦ lKijS ♦ II; Tj • OL 01 sec. Tj * Oi 02 sec 


O 1 PROCESSORS 
o 2 PROCESSORS 



vs; 


Figure 5. - Relative error (or second order system 
l/lrjS ♦ 1 Mt ? S ♦ ll t tj • ft 002 sec, • 0. 02 sec. 
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Si 
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O FIRST ORDER SYSTEM, 

T*aoe sec 

O SECOND ORDER SYSTEM, 

Ti • 0. 01 sec, Tj • tt 02 sec 
O SECOND CRIDER SYSTEM, 

t, - 11002 sec, t,- 0.02 sec 
O FOURTH ORDER SYSTEM, 

N T min -a0®fflsec 

o o 

\ 

. 41 - o \ 


\ 
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Figure 6. - Maximum stable time step with respect 
to minimum time constant of system for given 
number of processors. 
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Figure 7. - Effect d number of processor! on enptote mode* response to » *3 percent stop tot luto flow tor L 0 sec. 
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